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ABSTRACT 

We test the Secondary Infall Model (SIM) by direct comparison with the results of N- 
body simulations. Eight cluster-size and six galactic-size dark matter haloes have been 
selected at z = and re-simulated with high resolution. Based on their density pro- 
files at the initial redshift, we compute their evolution by the SIM, assuming a simple 
prescription for the angular momentum. A comparison of the density profiles obtained 
by the SIM and the numerical experiments at z = 5, 1 and shows that, for most 
of the haloes at most epochs, the SIM reproduces the simulated mater distribution 
with a typical fractional deviation of less than 40 per cent over more than six order of 
magnitudes in the density. It is also found that, within the SIM framework, most of 
the diversity in the shape of the density profiles at z = arises from the scatter in the 
primordial initial conditions rather than the scatter in the angular momentum distri- 
bution. A crude optimization shows that a similar degree of agreement is obtained for 
galactic and cluster haloes, but the former seem to require slightly higher amounts of 
angular momentum. Our main conclusion is that the SIM provides a viable dynamical 
model for predicting the structure and evolution of the density profile of dark matter 
haloes. 
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1 INTRODUCTION 

In the hierarchical clustering paradigm, cosmological struc- 
ture evolves from a primordial density perturbation field via 
gravitational instability. One of the basic outcomes of the 
process is the formation of virialized objects, usually referred 
to as dark matter (DM) haloes. The problem of the collapse 
and virialization of DM haloes can be easily formulated, as 
it depends on a single, relatively simple, physical process: 
the dissipationless gravitational interaction of a system com- 
posed of a large number of point-like particles. 

Yet, the long range and non-linear nature of the grav- 
itational N-body problem has made it extremely elusive 
to rigorous analytical treatment. On a phenomenologi- 
cal level, however, tremendous progress has been achieved 
by means of numerical experiments. After roughly three 
decades of N-bod y simulatio ns (starting with Aarse th et all 
[l979 : van Albada fc van Gor kom 1977; White 197a among 
others), this technique has reached the level at which differ- 
ent algorithms and numerical codes have converged to yield 
an overall consistent description of the formation, evolu- 
tion and internal structure of DM haloes jKnebe et alj200l 
I2OOJ) . 

A major breakthrough was the suggestion by 
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iNavarro et all Jl996l. Il997l) (hereafter NFW) that the den- 
sity profile of simulated DM haloes can be fitted by a simple 
analytical function. 



p(r) 



4ps 



(r/rs)(l + r/rs)2' 



(1) 



where the two free parameters ps and represent a char- 
acteristic density and radius of the halo. There is gen- 
eral consensus within the N-body community that the so- 
called NFW profile is able to provide a reasonably good fit 
to the numerical results in a variety of cold dark matter 
(CDM) cosmologies, although some doubts have been cast 
on the exact value of the logarit hmic slope at the centre 



(the so-called density cusp, see e.g. [Moore et al J|l99^. 
Ghigna et alJll998l |200^ iFukushi ge fc Makindll997l 
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2004; Na varro et alJ 12004^ . as well as on the degree of 
'universality' of the fit, i.e. its dependence on the un- 
derlying cosmological model, the mass accretion history 
of the halo, or it s environment (e.g. Ijing fc Sutc 
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200E 



iKlvDin et al.' '2001; Ricott] l2003l : lAvila-Reese et al 
Paulbetsch et al. 2006). 

Observationally, the DM density profile must be in- 
directly inferred. In ri ch clusters of galaxie s, estimates 
based on X-rays (e.g. IVoigt fc FabianlbOOel and refer- 
ences therein) or gravitational lensing (e.g. iDahle et alJ 
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I2OO3I: iGavazzi et al Bartelmann fc Meneehettil 12004 

iDalal fc Keeto J2003l) seem to indicate that the mass distri- 
bution is indeed well described by equation 1^3 . On galac- 
tic scales, observations are much more difficult to interpret. 
Dwarf spiral and low surface brightness galaxies are most 
likely dominated by DM in their centres, but the mass distri- 
bution inferred from rotation curves sugg ests a constant DM 
dens ity core rather than a steep cusp fe.g. lFlores fc Primackl 
Il994 : Moore. 1994) . Recent analyses show that observational 
data may actually be consistent with steeper profiles once 
the effects of inclination, non-circular orbits and triaxiality 
of the dark matter haloes are a ccounted for ||Havashi ct al. 
I2OO4I: iHavashi fc N avarrG"2006^ . bu t the controve rsy is still 
unresolved (e.g. ICentile et al. 20oi Ide Blokllioosll . 

It is therefore not clear whether the existence of a more 
or less 'universal' density profile is supported by observa- 
tional data, nor whether the possible disagreement is due to 
our lack of an understanding of the very complicated process 
of galaxy formation or it calls for an overall revision of the 
the standard model of cosmology. In any case, the qualita- 
tive shape of the density profile (shallower than isothermal 
in the inner parts and steeper in the outer), as well as (to 
some extent) its 'universality' constitute a very robust pre- 
diction of the currently accepted DM structure formation 
paradigm. 

Unfortunately, the NFW fitting formula provides a good 
phenomenological description of the density profile of simu- 
lated DM haloes, but it does not provide a physical under- 
standing of its origin. Ideally, one would have liked to have a 
self-consistent analytical model that enabled the calculation, 
from first principles only, of the density profile expected for 
a given cosmology, halo mass, formation time, environment, 
or whatever physically meaningful parameter that is found 
to play a relevant role. 

No such model exists as yet, but we are now close to 
having one. Much of the analytical work on t he formation of 
DM h aloes is based on the seminal paper of iGunn fc Gotd 
il972h on the dissipationless collapse of a spherical homo- 
geneous perturbation in an otherwise expanding Friedmann 
universe. This was followed by Gunn (1977), who consid- 
ered the collapse of a spherical inhomogeneous perturbation 
and the role of shell crossing. It also studied the case of 
the secondary infall, namely the late infall of shells onto an 
already collapsed and virialized perturbation, and showed 
that a self-similar solution can be found upon the use of 
adi abatic i n varian ce . 

iGunnI (|l973) was f ollowed by two differ e nt lin es of 
research. On one h and, iFillmore fc GoldreichI lll984) and 
iBertschingeJ 1^8^ independently found self-similar solu- 
tions of the collapse of scale-free spherical perturbations 
in an Einstein-de Sitter universe. On the other hand, 
iHoffman fc ShahamI 1^8^ analyzed the dependence of the 
structure of proto-haloes on the primordial power spectrum 
and the way it affects the density profile of virialized DM 
haloes. 

The basic predictions of the secondary infall model 
(SIM) were qualitatively confirmed by N-body simulations 
(e.g. Quiim et al. 1986: .Efstathio u et al. 1988; Crorie et alj 
Il994) . which prompted further study and extension of the 
SIM, mostly focusing on the nature of the initial condi- 
tions, on the improvement of the dynamical model and 
on its cosmological implications (e.g. iRvden fc Gunnlll987l : 



HoffmanI Il988l: iRvdenl Il988allbt IZaroubi fc HoffmanI [l993 : 
Lokas fc Hoffmanll2000l) . A common assumption made in 
these studies was that the halo particles are moving along 
radial orbits, and the generic result that emerges is that the 
inner density profile is roughly given by p oc . Thus, a 
cuspy density profile is a generic outcome of the SIM. Ac- 
tually, the predicted logarithmic slope is even steeper than 
the one measured in N-body experiments. 

The next major improvement of the SIM was the in- 
troduction of non-radial motions, namely a distribution of 
angular momentum of individual particles, not necessarily 
implying, on average, a total an gular momentum of the 
halo. A number of authors (e.g. [ White fc Zaritskvl Il992l : 
IR vdenlll99.?l: lAvila-R.eese et "Il] ll99?l " have pointed o ut that 
non-ra dial motions fiatten the inner density profile. iNusserl 
i200ll) extended the self-similar solutions to include a dis- 
tribution of angular momentum, and several recent studies 
show that, by introducing enough angular momentum, the 
density profil e predicted by the SIM can have a p oc 
density cusp llHiotelij 120021: iLe Delliou fc HenriksenI l2003l : 
lAscasibar et al.ll2004l: IWilliams et alJl2004l: iLu et al.l l2006h. 

Within the SIM framework, there are two key ingredi- 
ents that determine the shape of the virialized density pro- 
file: the initial density profile of the proto-halo and the an- 
gular momentum distribution. The primordial density has 
been determined either by employing the statistical prop- 
erties of Gaussian random field s (JIofKnan^.^.SIialMm.,J.98a; 
iBardeen et all Il986l: iHoffmanI [19881: lAscasibar et alJ 1200*11 
or by means of the mas s accretion history of DM haloes 
jAvila-Reese et alJ Il998l: iNusser fc ShethI Il999l: iLu et alJ 
l2Q0^)~The angular momentum distribution has always been 
considered to be a free parameter, d etermined so as t o yield 
best agreement with simulations. lAscasibar et al.l (j2003) 
pushed the study of the SIM one step further by making 
a detailed comparison between the density profile predicted 
by the model with the actual numerically simulated DM 
haloes rather than the NFW fit. These authors assumed the 
initial density profile to be given by the average expecta- 
tion around a local maximum of the primordial perturbation 
field, determined by the peak height and its smoothing scale. 
Optimizing over these two parameters and the angular mo- 
mentum distribution, the resulting density profiles were in 
close agreement with the simulated ones, with an accuracy 
comparable to that of the best-fitting NFW profile. 

The aim of the pr esent paper is to go one step beyond 
lAscasibar et al.l ll2004l) in order to assess the validity of the 
SIM as a tool to understand the physical origin of the DM 
density profile. More precisely, we attempt to recover the 
density profile of actual simulated haloes, but instead of as- 
suming a functional form for the initial conditions, we iden- 
tify the primordial peak in the initial snapshot of the sim- 
ulation and use the density profile around that point as an 
input for the SIM. We calculate the entire dynamical evo- 
lution of the halo upon assuming only one free parameter 
(the angular momentum distribution of the DM particles), 
and compare the density profile thus obtained with the full 
N-body simulation. We consider this as the ultimate test of 
the validity of the SIM, at least in the framework of the 
CDM-like cosmology. 

The paper is structured as follows: The basic principles 
of the SIM are reviewed in § [5| Numerical simulations are 
described in § [3 and the comparison of the SIM-calculated 
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density profiles with the results of the full N-body simula- 
tion is presented in § 0| A discussion and summary of the 
results is given in § |3 Technical details of the SIM and its 
implementation are given in Appendix IXI 



2 SECONDARY INFALL MODEL 

It is well known from cosmological N-body simulations that 
halo formation proceeds through a series of violent merger 
events involving smaller substructures. Indeed, major merg- 
ers seem to play an impo rtant role in shaping the st ructure of 
DM haloes (see e.g.lSv er & White 1998; Salvador-S ole et alJ 
Il998l: iManriaue et alJl2nn3. : .Romano-Diaz et alii2nO^ and 
therefore one would naively expect any model based on 
spherical symmetry to be hopelessly irrelevant and inade- 
quate. 

The SIM is not only based on the assumption of spher- 
ical symmetry. It further assumes that the complicated pro- 
cess of shell crossing, in which particles, represented by 
spherical shells, do not conserve their individual energies, 
admits a simple adiabatic invariant. This allows one to by- 
pass the need for a full self-consistent, one-dimensional, 
spherically symmetric, non-linear calculation of the dynam- 
ics of the collapsing shells, and use a semi-analytical method 
to calculate the equilibrium structure of DM haloes (for 
an interesting alternative a pproach, see the recent work by 
ISanchez -Conde etalJIgOOd) . Adiabatic invariance was intro- 
duced bv lGunnI Jl977l) .' who applied it to the case of self- 
similar collapse a nd this was extended to th e general non- 
power law case bv lZaroubi fc Hoflmanl Jl993l). Angular mo- 
mentum was first included by |^^^^^^^inn| I^S^) , and 
then in a more rigorous way bv lNusseJ (|200J^ . 

A very brief description of the SIM is given here and a 
detailed account of the SIM dynamical model and its imple- 
mentation is given in Appendix^ The physical model that 
underlies the SIM is based on an assumed primordial den- 
sity profile of the proto-object. Given that, the trajectory 
of a given mass shell is followed analytically until it reaches 
its turn-around radius. The shell does not experience shell 
crossing before it turns around, its energy is conserved and 
the calculation is exact. From that point, one needs to resort 
to the assumption that the shell trajectory admits an adi- 
abatic invariant, whereby the the maximum radius reached 
by the shell times the enclosed mass is approximately con- 
served. The transformation from the shell trajectories to a 
density profile is done by calculating the amount of time a 
given shell spends within a given radial interval. The angu- 
lar momentum of each shell is determined by a free param- 
eter, rj, that measures the square of the angular momentum 
of dark matter particles in terms of the maximal possible 
value and controls both the minimum radius reached by the 
shell and its time-radius relation. The physical basis for in- 
voking adiabatic invariance is that, as a given shell turns 
around and settles into equilibrium, it constitutes only a 
small perturbation to shells that have already collapsed and 
virialized. 

It follows that the evolution of the resulting DM halo 
is completely determined by its initial density profile and 
the intrinsic angular momentum distribution. In the ab- 
sence of a detailed information on the primordial struc- 
ture of a given halo, one must resort to statistical means. 



Two main approaches have been used here. One relies on 
the assumption that DM haloes are seeded by local den- 
sity maxima of the primordial perturbation field, which is 
assumed to be Gaussian. The mean density profile around 
a local density max imum can then be readily calculated 
jBardeen et alJll986D and use d to set the initial density pro- 
file of the proto-hal o (Hoff man &: Shahamlll985l : iHoffmaiil 
1198^1; lAscasibaL et alJI2004ll . Alternatively, one can use the 
mass accretion history of haloes, together with the spheri- 
cal top-hat model, to compute the mean density profile of 
proto-haloes llNusser fc Shethll999l : lAvila-Reese et al.ll99j : 
iLu et alJliooai ! 

The other key element of the SIM is the angular mo- 
mentum distribution of the DM particles. The amount of 
time a particle spends near the centre is controlled by its 
angular momentum, and thereby this quantity affects the 
'cross-talk' and en ergy exchange between t he sh ells. For pure 
radial orbits, the iFillmore fc GoldreichI (Il984l) self-similar 
solution is recovered almost independently of the initial con- 
dition s, and an de nsity cusp is obta ined both numeri- 
cally fHu ss et alJll99gri and analytically iLokas fc HoffmanI 
2000J . For the case of circul ar orbits, the turn around r adius 
density profile is recovered jHoffman fc Shaham|[l98il . 

The methodology of the present paper is to calculate 
the evolution of the density profile of individual DM haloes 
by the SIM and compare it with the simulated profiles. The 
initial conditions of the haloes are provided by the simu- 
lations themselves, rather than by some general considera- 
tions. Namely, the simulations are playing here a dual role, 
as they calculate the evolution of the selected objects, but 
also provide a mapping from the final virialized structures to 
their initial conditions. Here we improve on earlier studies of 
the SIM by testing the time evolution of the haloes and by 
replacing an assumed parametric fit to the initial conditions 
with the actual ones. 

As in previous work, we simply assume an angular mo- 
mentum distribution, and the empirical relation IIA14t has 
been invoked between the specific angular momentum of DM 
particles, the mass enclosed by the corresponding spherical 
shell and its turn-around radius. The parameter rj controls 
the amount of angular momentum injected, with the ex- 
treme values ri = and ri = 1 corresponding to purely ra- 
dial and circular orbits, respectively (see Appendix ^ for 
details). 



3 N-BODY SIMULATIONS 

We u se the Adaptive Refinement Tree code jKravtsov et alJ 
Il997tl to follow the evolution of structure in the concordance 
ACDM cosmology («,„ = 0.3, = 0.7, h = 0.7, as = 0.9). 

We are interested in the evolution of cluster-sized as 
well as galaxy-sized haloes. To construct suita ble initial con- 
dition s, we use the multiple-mass technique iKlvpin et al.l 
I2OOII) . In a first step we created an unconstrained random 
realization on an iV = 1024^ or iV = 2048^ grid. The ini- 
tial displacements and velocities of the particles were calcu- 
lated using all waves ranging from the fundamental mode 
k = 2tt/L to the Nyquist frequency kny = 2tt/L x iV^''^/2. 
Initial conditions at lower resolution were produced by merg- 
ing those small-mass particles to get in total 128^ particles 
and assigning to merged particles a velocity and a displace- 
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ment of one of the small-mass particles. The whole box of 
size was first simulated at this low resolution, 
starting at redshift zic = 50. 

Eight clusters have been selected from this low- 
resolution simulation, and the multiple-mass technique was 
used to set up high-resolution initial conditions. Namely, a 
Lagrangian region corresponding to a sphere of radius equal 
to two virial radii around each halo at z = was sampled 
with particles of mass rrip — 3.16 x 10* h~'^M0, correspond- 
ing to an effective number of 512^ particles in the box. The 
high mass-resolution region was surrounded by layers of par- 
ticles of increasing mass. The force resolution reached is 
2Ah~^'kpc (two times the size of the highest refinement 
level cell). 

To study the galaxy sized haloes, we have selected a 
filamentary region within a box of 80/i~^Mpc size. With 
the same resimulation technique as described above, this re- 
gion has been simulated with 150 million particles, which 
corresponds to an effective number of 2048"^ particles. Thus, 
the mass resolution is 5.0 x lO^h'^M©, and force resolu- 
tion reaches 0.6 kpc. Six galactic-scale haloes have been 
selected. 

For each object, we find the density maximum at 2: = 
by iteratively computing the centre of mass, starting with a 
reasonable initial guess and a search radius of 0.1 Mpc. 
Once convergence is reached, the search radius is reduced 
by ten percent, and the process is repeated until the sphere 
encloses less than 10^ particles. 

Then, the positions of the lO** particles closest to the 
density maximum at z = are traced back to the initial 
conditions, and we use only those positions to locate the 
primordial peak ai z — zic- It is interesting to note that, in 
many cases, the particles near the centre at 2 = belong to 
different objects at zic, which merged at some point during 
the history of the halo. Our procedure identifies the density 
peak corresponding to the most massive progenitor, which is 
a well defined entity, even at such a high redshift. Once the 
peak is located, we recompute its position, this time taking 
into account all particles at zic- The difference in position 
is usually not very large, but it has a significant impact on 
the shape of the primordial density profile near the centre. 

Finally, we trace the 10* particles closest to the density 
maximum at zic and locate the descendants at 2: = 5, 1 and 
0. We find the density maximum at those redshifts, taking 
into account all particles in the corresponding snapshots, 
and compute the simulated density profiles. 



4 RESULTS 

The dynamical evolution of the selected haloes has been cal- 
culated by both the SIM and the N-body simulation. By con- 
struction, the initial profiles coincide at zic, and we compare 
the results obtained by both methods at different redshifts 
(namely z — 5, 1 and 0). As explained in Section|5|and Ap- 
pendix^ angular momentum of the spherical shells in the 
SIM is determined by the 77 parameter. For each object, the 
SIM has been applied with a range of angular momentum 
values given by 77 = 0.10, 0.11, 0.12, 0.22. 

A first comparison of the SIM and N-body profiles is 
presented in Figure where all the SIM density profiles 
have been computed for 77 — 0.15. When this parameter 



is allowed to vary from object to object, a slightly better 
agreement can be obtained in the inner regions. However, 
it is not our aim to fit the N-body data, but to predict the 
evolution of the density profile, using only information at 
zic- Therefore, no attempt has been made here to make a 
quantitative determination of best-fitting rj for every object. 

A qualitative estimate may nevertheless shed some light 
on the most plausible values of ry, the associated scatter, 
and perhaps the dependence on mass, environment, or other 
factors. Figure 121 shows again the comparison between SIM 
and N-body profiles, but in this case the value of 77 has been 
individually optimized by visual inspection. In some cases, 
it might be possible to obtain an even closer correspondence 
between both methods if we also let rj vary with time, but 
in general terms, our results suggest that the introduction 
of such an additional degree of freedom is not necessary. 

It is interesting to compare here the galactic- and 
cluster-scale haloes. The size of the two samples is very small 
(six and eight objects) and no attempt is made here for any 
formal statistics, yet we find that in general galactic haloes 
are better fitted by higher values of rj. The G4 halo con- 
stitutes an enigmatic case. At z=0, the SIM provides an 
excellent fit to the simulated profile for an extremely high 
value of 77 = 0.60. A close inspection of the halo reveals that 
it is in the process of a major merger, and yet the agreement 
at 2 = is remarkably good. 

We plot in Figure |3 (for a fixed rj = 0.15) and Fig- 
ure 21 (for individually optimized values) the ratio between 
the enclosed mass profile predicted by the SIM and the one 
calculated by the N-body simulation. The mean and the 
standard deviation of MsiM/MNbody are shown in Figure |S] 
where the statistics is calculated over all haloes, galactic 
and cluster size, normalizing the radius by the turn-around 
radius of each object. The SIM calculations are done for 
a fixed 77 (upper row) and for optimized 77 per individual 
haloes (lower row). In both cases, the plots show very little 
bias (~ 20 per cent) and a scatter not larger than 40 per 
cent over most of the radial range. 

Finally, we are also interested in understanding the 
physical origin of the diversity of density profiles found in 
the N-body experiments. In terms of the SIM formalism, we 
would like to know if this diversity arises from scatter in 
the primordial conditions, in the specific angular momen- 
tum distribution of each halo, or both. The range of profiles 
exhibited by the simulated objects at 2: = is shown on 
the leftmost panel of Figure |S| The middle panel shows how 
the initial conditions are affecting the shape of the density 
profiles by grouping together all the SIM-calculated density 
profiles, evaluated at z = 0, for a fixed 77 = 0.15. The role of 
the angular momentum is illustrated on the rightmost panel, 
where the density profile predicted by the SIM for cluster 
C7 at 2 = is plotted for 77 = 0.10, 0.15 and 0.20. The cases 
?7 = (radial orbits, expected to yield an cusp), 77 = 1/2, 
and 77 = 1 (circular orbits, reflecting the turnaround density 
profile, and hence expected to result in a core structure) are 
given for reference. From Figure |5] we conclude that most 
of the diversity in the shapes of the virialized DM haloes is 
contributed by the variability in the initial conditions. The 
distribution of angular momentum, at least within the range 
of 77 = 0.15 ± 0.05 that seems to describe our sample, plays 
only a secondary role. 
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Figure 1. Cumulative overdensity profiles, A{r) 



3M(r) 



r [Mpc/h] 

1, of our six galactic haloes (designated Gl to G6) and eight clusters 



(CI to C8). Dashed and solid lines display the results of the SIM and the N-body experiments, respectively. From top to bottom, they 
correspond to 2 = 0, 1, 5, and 50 (the initial time of the simulations). The SIM is applied with the same angular momentum distribution 
for all haloes, given by r; = 0.15. The location of the virial (left) and turnaround (right) radii at the present epoch is indicated by the 
small vertical lines near the x-axis. 
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Figure 



r [Mpc/h] 

2. Same as Figure but with the optimal value of r} for each cluster, as indicated in the individual frames. 
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Figure 3. Ratio between tlie cumulative mass profiles, M(r), calculated by the SIM and the N-body experiments, plotted for z = 5 
(dotted lines), 1 (dashed lines) and (solid lines). The SIM is applied with the same angular momentum distribution for all haloes, given 
by r; = 0.15. The location of the virial (left) and turnaround (right) radii at the present epoch is indicated by the small vertical lines 
near the a;- axis. 
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4. Same as Figure^ but with the optimal value of r} for each cluster, as 



indicated in the individual frames. 
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Figure 5. Mean (solid line) and standard deviation (dashed lines) of the ratio between the enclosed mass predicted by the SIM and 
that computed by the N-body simulation, averaged over all (galactic- and cluster-size) haloes, scaled by their turn-around radius. Results 
obtained for fixed {rj = 0.15) and individually-optimized rj are plotted in the upper and lower panels, respectively. The three columns 
correspond to redshifts 2 = 0, 1 and 5. Horizontal dotted lines indicate a fra<;tional deviation of ±40 per cent. 
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Figure 6. Diversity in the shape of the density profile. Left ijaiicl: radial density profile, scaled by the turn-around radius, of all 
simulated objects at the present epoch. Middle panel: profiles calculated by the SIM, assuming a fixed angular momentum distribution 
given by r; = 0.15. Right panel: Effect of ajigular momentum in the SIM. Talcing object C7 as reference, the solid (?? = 0.15) and two 
dashed lines (77 = 0.10,0.20) encompass the typical values of rj found in our study. The extreme cases r] = 0, 0.5, and 1.0 (dotted lines) 
and power laws and (small solid lines) are also shown for reference. The comparison shows that, within a typical range of r], the 
diversity in the shapes of the density profile at 2 = is dominated by variations in the primordial structure rather than by the amount 
of angular momentum of each halo. 
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5 DISCUSSION 

The main result of the present work is that the secondary 
infall model provides a valid theoretical framework for cal- 
culating the structure and evolution of DM haloes in an 
expanding universe, and that its predictions with respect to 
the density profiles are in close agreement with full N-body 
simulations. Comparing the SIM and simulated cumulative 
density profiles over more than six orders of magnitude, we 
find the typical discrepancy to be better than 40 per cent. 
This level of agreement extends to the time evolution of the 
density profile and it is not limited to the final time snap- 
shot only. Within the SIM framework, most of the diversity 
in the density profiles of DM haloes is contributed by the 
scatter in the primordial profiles rather than the scatter in 
the angular momentum. 

Galactic-size haloes are less isolated objects than cluster 
haloes, in the sense that they have more similar-mass com- 
panions. The simulated galaxies have been selected from a 
typical filamentary region and naive reasoning would lead 
one to expect these objects to be less suitable for the appli- 
cation of the SIM compared with clusters. This is not the 
case. The agreement between the N-body and SIM density 
profile of galactic haloes is roughly as good as for clusters. 
The only difference is that the crude optimization over the 
angular momentum parameter, 77, finds that galaxies require 
slightly higher values of 77, consistent with the fact that the 
galactic haloes are more torqued and hence should have ac- 
quired more angular momentum. This in turn gives further 
support to the SIM in the sense that the value of rj is not 
a mere free parameter but is related to the actual angular 
momentum of haloes. 

The galactic-scale halo G4 deserves a special attention. 
The SIM-predicted profile with = 0.15 provides a very bad 
fit to the simulated density profile at z — 0. Yet, for the very 
high value of 77 = 0.60, a very good fit is obtained for the 
present epoch. Such a high value constitutes a remarkable 
exception, and no other object has a best-fitting rj anywhere 
close to 0.60. Visual inspection of the G4 halo reveals a 
major merger in progress at 2: = 0, and the halo is clearly 
not in virial equilibrium. In spite of its dynamical state, the 
density profile of the halo is extremely well fitted by the 
SIM. Is such an agreement a mere coincidence, or does the 
SIM actually reproduce the internal structure the system? 
This intriguing question has led us to further investigate 
the structure of the G4 halo, and we found that it actually 
displays an extremely high amount of angular momentum 
due to the merger event. 

Although the present work focuses on the density pro- 
file, it would be desirable to extend the formalism outlined 
in Appendix^so as to compute their full dynamical struc- 
ture of DM haloes; more precisely, the radial and tangen- 
tial components of the velocity dispersion. Derived quan- 
tities, such as the angular momentum distribution or the 
anisotropy profile, can be readily calculated, and their de- 
pendence on the assumed value of ry would be extremely 
helpful in validating our simple prescription for angular mo- 
mentum quantitatively. 

Moreover, one could use the upgraded formalism to 
address additional questions raised by numerical N-body 
simulations, such as the relation between the logarithmic 
slope of the density profile and the anisotropy of the ve- 



locity dispersion jHansen fc Moor3l2006h . or the power-law 
radial dependence of the coarse-grained phase-space density 
('Taylor fc Navarro! I2OO1I lAscasibar et al.ll2004l: iRasia et all 

2004). This quantity is somewhat analogous to the 'entropy' 
of the DM, formally defined so as to mimic the ideal gas 
thermal entropy, and it has indeed been found that they 
are closely related in ad iabatic gasdynamical simulations 
jFaltenbacher et alJ '2006\ 

Eventually, the final goal would be to achieve a totally 
self-consistent prescription to evaluate the angular momen- 
tum distribution of DM haloes, which of course can only 
be achieved through an understanding of its physical ori- 
gin. Although strict spherical symmetry would imply that 
the angular momentum of individual particles must be ex- 
actly conserved (and hence it should be imprinted by the 
initial conditions), several processes may be responsible for 
the conversion of radial into tangential motions in a more 
realistic case. 

For instance, analytical and numerical studies have 
shown that a predominantly radial collapse is unstable to 
the growth of tangential motions, resulting in the so-called 
radial orbit instability (see e.g. Huss et al. 1 9991, and ref- 
erenc es therein), and it has been conjectured llBarnes et alJ 
l2006h that this mechanism could actually set the length scale 
Vs of DM haloes, which marks the transition region from al- 
most isotropic orbits to radial ones. 

A seemingly very different point of view is provided by 
cosmological simulations, which are assumed to incorporate 
all the physics relevant for the formation of DM haloes. 
Haloes in numerical experiments evolve through phases of 
violent 'major merger' events in which objects of compara- 
ble masses merge to form a yet bigger halo. A detailed anal- 
ysis shows that the NFW-like structure, and in particular 
the angular momentum distribution within rs, is exclusively 
determined by major mergers, and is hardly affected by 
the quiescent phases that follow iSalvador-Sole et "aPbOOSl : 
iRomano-Diaz et alJl20o3 . and references therein). 

To summarize, two very different frames of reference for 
DM halo formation end up predicting very similar density 
profiles. Assuming that this is not a coincidence, one should 
look for an underlying physical process that unites these 
two different approaches. Our results suggest that angular 
momentum might play an important role in that respect. 

Yet, the present study points to another key element 
and this is the initial configuration of the haloes. Figure |S| 
shows the diversity in the density profile of the virialized 
objects caused by the variation in their initial density struc- 
ture. It follows that the shape of the halo density profile 
cannot be discussed without an explicit reference to the na- 
ture of the initial conditions; a link that any complete theory 
of the collapse of DM haloes should incorporate. Although a 
full understanding is yet to be achieved, the study of the col- 
lapse and virialization of cosmic structures that commenced 
with iGunn fc Gotd lll973) has certainly provided us with 
many clues and a deep insight of the most relevant ingredi- 
ents behind the process. 

As a final note, we would like to suggest a potential 
practical application of the SIM as a complement (per- 
haps even substitute) of the multiple-mass technique. High- 
resolution re-simulations of individual DM haloes identified 
in low-resolution cosmological experiments require a signif- 
icant amount of CPU resources, and they can be practi- 
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cally performed only for a small subset of all objects. We 
have shown, however, that the SIM provides a valid frame- 
work for calculating the dynamical evolution of DM haloes, 
and it takes about 10 s to compute the density profile of a 
given object at a given epoch on a desktop PC. This gives 
rise to the interesting possibility of using the SIM to in- 
sert highly-resolved spherical haloes into large-scale, low- 
resolution cosmological simulations. In practice, one can use 
the low-resolution experiment to set up the initial conditions 
for all the objects and use the SIM to calculate their dynam- 
ical evolution, providing an efficient and relatively accurate 
way of generating extended catalogs of dynamically resolved 
clusters at a reasonable accuracy. 
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APPENDIX A: IMPLEMENTATION OF THE 
SPHERICAL INFALL MODEL 

Al Physics 

We consider the evolution of a spherical di stribution M(r), 
follow ing the prescription described in lAscasibar et alJ 
i2004l) . We use the initial comoving radius a; as a Lagrangian 
coordinate identifying a shell of matter enclosing a mass 



Before shell-crossing. 



47!" 3 
Mj: = —Q,m,pcX 



(Al) 



where Q.m denotes the current dark matter density and pc 
is the critical density. The physical location of the shell can 
be written in terms of the variable a{t), 



r{t) = xaia{t) 



(A2) 



where at = a{ti) is the expansion factor of the universe at 
an arbitrary initial time, ti. 

In a homogeneous universe, a{t) would simply be the 
cosmic expansion factor, i.e. a{t) = a{t)/ai Va;. A spherical 
perturbation with overdensity 



1 < 1 



(A3) 



can be obtained by slightly displacing the shells. To first 
order in Ai, 



r{ti) = xaia{ti) « xui 
while the velocity would be given by 



A. 



(A4) 



f{ti) = xaia{ti) « Hixai ( 1 ^ ) (A5) 



2Ai 



(A6) 



is a conserved quantity, analogous to the Newtonian specific 
energy of each shell. With our initial conditions, ei{x) ~ 
—5/6{Hixai)^Ai, yielding the equation of motion 



(A7) 



where Q.i = Q,mOL^ "^Ho/Hf and A^ = QAHo/Hf are the 
matter and vacuum energy densities at time ti, and 

K,^ -l + n,+A, + ^{4. + n,-2A.)^^. (A8) 

In fact, we typically obtain values of Ai ~ 0.3 in our 
study. Keeping only linear terms in A; results in errors of 
~ 25 per cent a.t t = ti, and we expect the m to grow as the 
system evolves. We will thus depart from lAscasibar et al.l 
(.2004^ and take the more accurate initial conditions 



a{t,) = (1 + AO'/' 
and compute Ki by solving 

Jo -JVLiOi-^ + A,q2 _ Xi Jo VTUo^+Jho? 

(AlO) 

that is, by imposing the correct initial time for every shell 
^ . The initial velocity can be easily obtained as 



(A9) 



da 



d(tO = Y/n,(l + A0i/3+A,(l-f A0-2/«-A-,. (All) 

For positive overdensities, each shell will expand slower 
than the average universe, and, for high enough Ai, will 
eventually turn around at some point _Rta = ^(?ta) and 
start collapsing, crossing the inner shells on its way towards 
the centre. After shell- crossing (which we approximate as 
immediately after turn-around), the appropriate integral of 
motion is 



£2(3;:) 



G 



r 2r 



(A12) 



where (j){r, t) is the gravitational potential of the halo and 
j is the average modulus of the specific angular momentum 
of the particles constituting the shell. We approximate the 
local potential by a power-law mass distribution 



GM(x) , 
[r] = I ^—'-dx : 



.7-1 



GM, / X y-^ GM^ 



(A13) 

where 7? is the maximum radius of the orbit, 7 is the local 
logarithmic slope of the mass profile, f{x) = ln(2;) for 7=1 
1) otherwise. Concerning angular 



and f{x) = 
momentum, we adopt the prescription 



(A14) 



where 77 is a free parameter from (radial) to 1 (circular 
orbits). For test particles orbiting a point mass (7 = 0), 
the eccentricity of the orbits would he e = 1 — r], while the 
pericentric radius would be fmi-a/R = (1 — e)/(l -\- e) = 



where Hi = ^ ^ma^'^ + -I- (1 — fim — Q,A)a'l'^ corre- 
sponds to the Hubble constant aX t = ti. In what follows, 
time will always be expressed in units of . 



^ For f2i ~ 1 and ~ 0, a good approximation at 

high redshift, the implicit equation lAlOi would reduce to 



3/2 (.Xi(l+Ai)-l/3 

Jo 



dx - 
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7)/(2 — 77). The equation of motion, written in terms of A ; 
r/R = a/aR, is 



Initially, R = i?ta, but the potential (ji{r) is not static. 
Outer shells that collapse later will add a certain amount 
of mass, Ms_dd{r), and the apocentre R will slowly move in- 
wards. If (/!>(r) changes slowly compared to the orbital period, 
both the radial action 



Jr ^ <j) rdr oc ^RM{R) 



(A16) 



and the angular action (the specific angular momentum j) 
should be conserved. We find the final apocentric radius, R' , 
by solving the implicit equation 



R! 



R M^ + Mm{R') 



(A17) 



numerically. Constant angular momentum implies -q' = 77, 
while the added mass should be taken into account by 



^^R_ 
Hi R'' 



(A18) 



In order to compute A/add (»■), we assume instantaneous 
phase mixing. After turn-around, each spherical shell turns 
into a density distribution with cumulative mass propor- 
tional to the fraction of time its constituent particles spend 
within r, 



dM,dd(r) = 



dM^ 
T 



"■/^dA 
A 



(A19) 



where dMx is the mass of the shell and T — ^^AX/\ its 
orbital period. Neglecting A^, 



Madd(r) ^ / r 



with 



/-[f(l-A-^)-/(A)]-^/^dA 
/o^[2(l-A-2)-/(A)]-i/2dA' 



(A20) 



(A21) 



A2 Numerical scheme 

First of all, the program computes the initial conditions from 
a given mass profile, M(r), obtained from a simulation snap- 
shot at redshift zic'- 



3M 



1 (1 + zic)ai 



(A22) 



In our case, the snapshot corresponds to zic ~ 50. The 
initial time for the SIM integration is set by the value of a^, 
which we choose, for the sake of simplicity, to coincide with 
the time of the simulation snapshot, i.e. ai = 1/(1 + zic)- 

Then, the evolution of logarithmically-spaced shells 
is integrated numerically, starting by the outermost, less 
overdense one. For shells still expanding at the final epoch, 
Oo, one simply computes ag — a{to) from 1A7II . The ensuing 
profile is, not surprisingly. 



The situation becomes more complicated after turn- 
around, where 



M{R') = M, + A'Ui{R') 



(A24) 



M{xaiao) = 



(A23) 



For these shells, we compute Rta,{Mx) ~ xaiata. by find- 
ing the zero of equation IIA7l l. Then we compute R'{Mx) 
according to llA17ll and add the contribution of this shell 
to Madd('') ^ using equation HAIQII . In order to obtain a 
smooth transition between the two regimes, the contribu- 
tion to Madd('^) of shells with turn-around time Tta > to/2 
is multiplied by a factor to /Tta. — 1. 

The process is repeated iteratively, decreasing the shell 
mass Mx by 1 per cent, until the initial radius is smaller 
than the innermost data point in the initial conditions. 
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